{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Baysian Excess Variance (Bexvar)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Bayesian Excess Variance (bexvar) is a statistical measurement of variability in Poisson-distributed light curves. Bexvar is a Bayesian formulation of excess variance. A brief summary of theoretical understanding of bexvar is given at the end of this tutorial. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `bexvar()` method implemented in Stingray, provides posterior samples of bexvar given a light curve data as input parameters. \n", "This tutorial is intended to give a demonstration of How to use `bexvar()` method implemented in Stingray.\n", "The method takes following input parameters. (Given here for completeness)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "  ```time``` : iterable, `:class:numpy.array` or `:class:List` of floats, optional, default ``None`` \n", "     A list or array of time stamps for a light curve. \n", "  `time_del` : iterable, `:class:numpy.array` or `:class:List` of floats \n", "    A list or array of time intervals for each bin of light curve. \n", "  `src_counts` : iterable, `:class:numpy.array` or `:class:List` of floats \n", "    A list or array of counts observed from source region in each bin. \n", "  `bg_counts` : iterable, `:class:numpy.array` or `:class:List` of floats, optional, default ``None`` \n", "    A list or array of counts observed from background region in each bin. If ``None`` \n", "    we assume it as a numpy array of zeros, of length equal to length of ``src_counts``. \n", "  `bg_ratio` : iterable, `:class:numpy.array` or `:class:List` of floats, optional, default ``None`` \n", "    A list or array of source region area to background region area ratio in each bin. \n", "    If ``None`` we assume it as a numpy array of ones, of length equal to the length of \n", "    ``src_counts``. \n", "  `frac_exp` : iterable, `:class:numpy.array` or `:class:List` of floats, optional, default ``None`` \n", "    A list or array of fractional exposers in each bin. If ``None`` we assume it as \n", "    a numpy array of ones, of length equal to length of ``src_counts``. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us start by importing the bexvar module" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "from stingray import bexvar" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now consider an example dataset." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "time = np.arange(0,8)*100\n", "counts= np.array([106, 87, 115, 148, 43, 129, 204, 87])\n", "time_del = np.ones(np.size(time))*100\n", "bg_counts = np.array([722, 696, 701, 721, 722, 703, 722, 695])\n", "bg_ratio = np.array([0.01474, 0.01158, 0.01214, 0.01308, 0.010877, 0.01177, 0.01058, 0.01138])\n", "frac_exp = np.array([0.37416, 0.21713, 0.37937, 0.50140, 0.11617, 0.39221, 0.64275, 0.31160])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Call bexvar function to get posterior distribution of bexvar." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `bexvar()` method uses [UltraNest](https://johannesbuchner.github.io/UltraNest/) python package to obtain the posteriors of bexvar. Ultranest gives a brief summary of log evidence (log(z)) and its uncertainties, and the parameter constraints. \n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "preparing time bin posteriors...\n", "running bexvar...\n", "[ultranest] Sampling 400 live points from prior ...\n", "[ultranest] Explored until L=-2e+01 [-20.4040..-20.4040]*| it/evals=3622/5046 eff=77.9595% N=400 \n", "[ultranest] Likelihood function evaluations: 5051\n", "[ultranest] logZ = -24.86 +- 0.0784\n", "[ultranest] Effective samples strategy satisfied (ESS = 1590.2, need >400)\n", "[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.47+-0.06 nat, need <0.50 nat)\n", "[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.08, need <0.5)\n", "[ultranest] logZ error budget: single: 0.09 bs:0.08 tail:0.01 total:0.08 required:<0.50\n", "[ultranest] done iterating.\n", "\n", "logZ = -24.856 +- 0.156\n", " single instance: logZ = -24.856 +- 0.093\n", " bootstrapped : logZ = -24.856 +- 0.156\n", " tail : logZ = +- 0.010\n", "insert order U test : converged: True correlation: inf iterations\n", "\n", " logmean : 0.350 │ ▁ ▁ ▁▁▁▁▁▁▁▁▂▃▄▅▆▇▇▇▆▅▄▃▂▁▁▁▁▁▁▁ ▁ ▁▁ │0.575 0.461 +- 0.020\n", " logsigma : 0.010 │▇▅▄▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁ ▁▁▁▁ ▁ ▁ │0.227 0.028 +- 0.018\n", "\n", "running bexvar... done\n" ] } ], "source": [ "\n", " bexvar_distribution = bexvar.bexvar(time=time, src_counts=counts, time_del=time_del, frac_exp=frac_exp,\n", " bg_counts=bg_counts, bg_ratio=bg_ratio)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can then plot the samples to visualize the posterior distribution of bexvar." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEICAYAAACuxNj9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY3ElEQVR4nO3debRlZXnn8e/PQnACBbkaoMACG0w7FlpBE2MWOCLYgENr0VFwaHHATmxcrTh0HBJ6kcGJNo1dKiJGQZSwRIU2SAftdItQYFkMghSDocoSCmxxDA349B/nvXIo7r371K17huJ+P2vtdfZ59rv3eWqfA8/d+9373akqJEmaywPGnYAkafJZLCRJnSwWkqROFgtJUieLhSSpk8VCktRpu3EnMCy77rprLVu2bNxpSNI249JLL721qqZmWna/LRbLli1j9erV405DkrYZSX442zJPQ0mSOlksJEmdLBaSpE4WC0lSJ4uFJKmTxUKS1MliIUnqZLGQJHW6396UtzWWHf+1ea9744mHLmAmkjQZPLKQJHUaWrFIckqSW5Jc0Rf7QpI1bboxyZoWX5bk133LPt63ztOSXJ5kXZKTkmRYOUuSZjbM01CnAh8DTpsOVNUrpueTfBC4va/9dVW1fIbtnAy8HvgOcC5wMHDewqcrSZrN0I4squpbwE9mWtaODl4OnD7XNpLsBuxUVRdVVdErPEcscKqSpA7j6rN4FnBzVV3bF9s7yXeTfDPJs1psD2B9X5v1LTajJMckWZ1k9aZNmxY+a0lapMZVLI7k3kcVG4G9qmp/4Djg80l22tKNVtWqqlpRVSumpmYckl2SNA8jv3Q2yXbAS4CnTceq6g7gjjZ/aZLrgP2ADcDSvtWXtpgkaYTGcWTxXODqqvrt6aUkU0mWtPl9gH2B66tqI/CzJM9o/RxHAV8eQ86StKgN89LZ04FvA49Lsj7J69qildy3Y/uPgLXtUtovAW+squnO8TcDnwTWAdfhlVCSNHJDOw1VVUfOEn/1DLGzgLNmab8aeOKCJidJ2iLewS1J6mSxkCR1slhIkjpZLCRJnSwWkqROFgtJUieLhSSpk8VCktTJYiFJ6mSxkCR1slhIkjpZLCRJnSwWkqROFgtJUieLhSSpk8VCktTJYiFJ6mSxkCR1slhIkjpZLCRJnYZWLJKckuSWJFf0xd6XZEOSNW06pG/ZO5OsS3JNkhf0xQ9usXVJjh9WvpKk2Q3zyOJU4OAZ4h+uquVtOhcgyeOBlcAT2jr/LcmSJEuAvwVeCDweOLK1lSSN0HbD2nBVfSvJsgGbHw6cUVV3ADckWQcc0Jatq6rrAZKc0dpetdD5SpJmN44+i7ckWdtOU+3cYnsAN/W1Wd9is8VnlOSYJKuTrN60adNC5y1Ji9aoi8XJwGOB5cBG4IMLufGqWlVVK6pqxdTU1EJuWpIWtaGdhppJVd08PZ/kE8BX29sNwJ59TZe2GHPEJUkjMtIjiyS79b19MTB9pdQ5wMokOyTZG9gXuBi4BNg3yd5JtqfXCX7OKHOWJA3xyCLJ6cCBwK5J1gPvBQ5Mshwo4EbgDQBVdWWSM+l1XN8FHFtVd7ftvAX4OrAEOKWqrhxWzpKkmQ3zaqgjZwh/ao72JwAnzBA/Fzh3AVOTJG0h7+CWJHWyWEiSOlksJEmdLBaSpE4WC0lSJ4uFJKmTxUKS1MliIUnqZLGQJHWyWEiSOlksJEmdLBaSpE4WC0lSJ4uFJKmTxUKS1MliIUnqZLGQJHWyWEiSOlksJEmdLBaSpE5DKxZJTklyS5Ir+mJ/neTqJGuTnJ3kES2+LMmvk6xp08f71nlaksuTrEtyUpIMK2dJ0syGeWRxKnDwZrHzgSdW1ZOBHwDv7Ft2XVUtb9Mb++InA68H9m3T5tuUJA3Z0IpFVX0L+MlmsX+oqrva24uApXNtI8luwE5VdVFVFXAacMQQ0pUkzWGcfRavBc7re793ku8m+WaSZ7XYHsD6vjbrW2xGSY5JsjrJ6k2bNi18xpK0SI2lWCR5N3AX8LkW2gjsVVX7A8cBn0+y05Zut6pWVdWKqloxNTW1cAlL0iK33ag/MMmrgRcBz2mnlqiqO4A72vylSa4D9gM2cO9TVUtbTJI0QiM9skhyMPB24LCq+lVffCrJkja/D72O7OuraiPwsyTPaFdBHQV8eZQ5S5K2sFgk2TnJkwdsezrwbeBxSdYneR3wMWBH4PzNLpH9I2BtkjXAl4A3VtV05/ibgU8C64DruHc/hyRpBDpPQyW5EDistb0UuCXJ/66q4+Zar6qOnCH8qVnangWcNcuy1cATu/KUJA3PIEcWD6+qnwEvAU6rqqcDzx1uWpKkSTJIsdiu3e/wcuCrQ85HkjSBBikWHwC+Tu8O60taB/S1w01LkjRJOvssquqLwBf73l8PvHSYSUmSJkvnkUWS/ZJcMD0gYJInJ3nP8FOTJE2KQU5DfYLegH93AlTVWmDlMJOSJE2WQYrFQ6rq4s1id83YUpJ0vzRIsbg1yWOBAkjyMnpjOUmSFolBxoY6FlgF/G6SDcANwCuHmpUkaaIMcjXU9cBzkzwUeEBV/Xz4aUmSJsmsxSLJjMN5TD/VtKo+NKScJEkTZq4jix1HloUkaaLNWiyq6v2jTESSNLkGuSlvnyRfSbIpyS1JvtyG/JAkLRKDXDr7eeBMYDdgd3pDf5w+zKQkSZNl0JvyPltVd7Xp74AHDTsxSdLkGOQ+i/OSHA+cQe/GvFcA5ybZBaDviXaSpPupQYrFy9vrGzaLr6RXPOy/kKT7uUFuytt7FIlIkibXIFdDLUlyWJI/SXLc9DTIxpOc0q6guqIvtkuS85Nc2153bvEkOSnJuiRrkzy1b52jW/trkxw9n3+oJGn+Bung/grwauCR9G7Um54GcSpw8Gax44ELqmpf4IL2HuCFwL5tOgY4GXrFBXgv8HTgAOC90wVGkjQag/RZLK2qJ89n41X1rSTLNgsfDhzY5j8DXAi8o8VPq6oCLkryiPbs7wOB86c70pOcT68AefmuJI3IIEcW5yV5/gJ+5qOranqI8x8Dj27zewA39bVb32KzxSVJIzLIkcVFwNlJHkDvaXkBqqp22toPr6pKUlu7nWlJjqF3Cou99tproTYrSYveIEcWHwJ+n97NeTtV1Y5bWShubqeXaK+3tPgGYM++dktbbLb4fVTVqqpaUVUrpqamtiJFSVK/QYrFTcAVrS9hIZwDTF/RdDTw5b74Ue2qqGcAt7fTVV8Hnp9k59ax/fwWkySNyCCnoa4HLkxyHnDHdHCQ51kkOZ1eB/WuSdbTu6rpRODMJK8Dfsg9N/2dCxwCrAN+Bbymfc5Pkvw5cElr9wHvGpek0RqkWNzQpu3bNLCqOnKWRc+ZoW3Re4TrTNs5BThlSz5bkrRwBrmD2+dabIFlx39t3uveeOKhC5iJJC2czmKRZAp4O/AE+kabrapnDzEvSdIEGaSD+3PA1cDewPuBG7mn/0CStAgMUiweWVWfAu6sqm9W1WsBjyokaREZpIP7zva6McmhwI+AXYaXkiRp0gxSLP4iycOBtwH/FdgJ+I9DzUqSNFEGuRrqq232duCg4aYjSZpEgzzP4q+S7JTkgUkuSLIpyStHkZwkaTIM0sH9/Kr6GfAieldC/SvgPw0zKUnSZBmkWEyfqjoU+GJV3T7EfCRJE2iQDu6vJrka+DXwpnaT3r8MNy1J0iTpPLKoquOBPwBWVNWd9Ab5O3zYiUmSJscgRxb0j/JaVb8Efjm0jCRJE2eQPgtJ0iI3a7FI8sz2usPo0pEkTaK5jixOaq/fHkUikqTJNVefxZ1JVgF7JDlp84VV9SfDS0uSNEnmKhYvAp4LvAC4dDTpSJIm0azFoqpuBc5I8v2q+t4Ic5IkTZhBroa6LcnZSW5p01lJlg49M0nSxBikWHwaOAfYvU1fabF5SfK4JGv6pp8leWuS9yXZ0Bc/pG+ddyZZl+SaJC+Y72dLkuZnkJvyHlVV/cXh1CRvne8HVtU1wHKAJEuADcDZwGuAD1fV3/S3T/J4YCW9Z4DvDnwjyX5Vdfd8c5AkbZlBjixuTfLKJEva9ErgtgX6/OcA11XVD+doczhwRlXdUVU3AOuAAxbo8yVJAxikWLwWeDnwY2Aj8DJ6RwELYSVwet/7tyRZm+SUJDu32B7ATX1t1reYJGlEBhlI8IdVdVhVTVXVo6rqiKr656394CTbA4cBX2yhk4HH0jtFtRH44Dy2eUyS1UlWb9q0aWtTlCQ14xwb6oXAZVV1M0BV3VxVd1fVb4BPcM+ppg3Ann3rLW2x+6iqVVW1oqpWTE1NDTF1SVpcxlksjqTvFFSS3fqWvRi4os2fA6xMskOSvYF9gYtHlqUkabAhyhdakocCzwPe0Bf+qyTLgaL3+NY3AFTVlUnOBK4C7gKO9UooSRqtzmKR5D1V9RdtfoequmNrP7Q9E+ORm8VeNUf7E4ATtvZzJUnzM9cQ5e9I8vv0rn6a5gi0krQIzXVkcTXwb4F9kvyv9v6RSR7XbqyTJC0Sc3Vw/xR4F72b4A4EPtrixyf5P8NNS5I0SeY6sngB8Gf07n34ELAW+GVVLdQNeZKkbcSsRxZV9a6qeg69K5M+CywBppL8U5KvjCg/SdIEGOTS2a9X1WpgdZI3VdUfJtl12IlJkibHIMN9vL3v7atb7NZhJSRJmjxbdAe3T8yTpMVpnMN9SJK2ERYLSVIni4UkqZPFQpLUyWIhSepksZAkdbJYSJI6WSwkSZ0sFpKkTmN5rKpmtuz4r8173RtPPHQBM5Gke/PIQpLUyWIhSeo0tmKR5MYklydZk2R1i+2S5Pwk17bXnVs8SU5Ksi7J2iRPHVfekrQYjfvI4qCqWl5VK9r744ELqmpf4IL2HuCFwL5tOgY4eeSZStIiNu5isbnDgc+0+c8AR/TFT6uei4BHJNltDPlJ0qI0zmJRwD8kuTTJMS326Kra2OZ/DDy6ze8B3NS37voWu5ckxyRZnWT1pk2bhpW3JC0647x09g+rakOSRwHnJ7m6f2FVVZLakg1W1SpgFcCKFSu2aF1J0uzGdmRRVRva6y3A2cABwM3Tp5fa6y2t+QZgz77Vl7aYJGkExlIskjw0yY7T88DzgSuAc4CjW7OjgS+3+XOAo9pVUc8Abu87XSVJGrJxnYZ6NHB2kukcPl9V/yPJJcCZSV4H/BB4eWt/LnAIsA74FfCa0acsSYvXWIpFVV0PPGWG+G3Ac2aIF3DsCFKTJM1g0i6dlSRNIIuFJKmTxUKS1MliIUnqZLGQJHWyWEiSOlksJEmdLBaSpE4WC0lSJ4uFJKnTOIco1wJadvzXtmr9G088dIEykXR/5JGFJKmTxUKS1MliIUnqZLGQJHWyWEiSOlksJEmdLBaSpE4WC0lSJ4uFJKnTyItFkj2T/GOSq5JcmeRPW/x9STYkWdOmQ/rWeWeSdUmuSfKCUecsSYvdOIb7uAt4W1VdlmRH4NIk57dlH66qv+lvnOTxwErgCcDuwDeS7FdVd480a0laxEZ+ZFFVG6vqsjb/c+D7wB5zrHI4cEZV3VFVNwDrgAOGn6kkadpY+yySLAP2B77TQm9JsjbJKUl2brE9gJv6VlvPLMUlyTFJVidZvWnTpmGlLUmLztiKRZKHAWcBb62qnwEnA48FlgMbgQ9u6TaralVVraiqFVNTUwuZriQtamMpFkkeSK9QfK6q/h6gqm6uqrur6jfAJ7jnVNMGYM++1Ze2mCRpRMZxNVSATwHfr6oP9cV362v2YuCKNn8OsDLJDkn2BvYFLh5VvpKk8VwN9UzgVcDlSda02LuAI5MsBwq4EXgDQFVdmeRM4Cp6V1Id65VQkjRaIy8WVfVPQGZYdO4c65wAnDC0pCRJc/IObklSJ4uFJKnTOPosNIGWHf+1ea9744mHLmAmkiaRRxaSpE4WC0lSJ4uFJKmTxUKS1MliIUnq5NVQ2mpeSSXd/3lkIUnqZLGQJHWyWEiSOlksJEmd7ODWWNk5Lm0bPLKQJHWyWEiSOlksJEmdLBaSpE4WC0lSJ6+G0jZra66kAq+mkrbENnNkkeTgJNckWZfk+HHnI0mLyTZxZJFkCfC3wPOA9cAlSc6pqqvGm5m2Zd7jIQ1umygWwAHAuqq6HiDJGcDhgMVCY7G1p8DmyyKlcdlWisUewE1979cDT9+8UZJjgGPa218kuWYEuW2pXYFbx53ELCY1t0nNC0acW/5yi5pP6n6b1LzA3B4z24JtpVgMpKpWAavGncdckqyuqhXjzmMmk5rbpOYF5jYfk5oXmNtctpUO7g3Ann3vl7aYJGkEtpVicQmwb5K9k2wPrATOGXNOkrRobBOnoarqriRvAb4OLAFOqaorx5zWfE3yabJJzW1S8wJzm49JzQvMbVapqnF+viRpG7CtnIaSJI2RxUKS1MlisRW6hiBJskOSL7Tl30myrMWfl+TSJJe312f3rXNh2+aaNj1qxLktS/Lrvs//eN86T2s5r0tyUpKMOLc/7strTZLfJFnelo1qv/1RksuS3JXkZZstOzrJtW06ui++1fttvnklWZ7k20muTLI2ySv6lp2a5Ia+fbZ8S/Pamtzasrv7Pv+cvvje7btf134L248qryQHbfY7+5ckR7Rlo9pnxyW5qn1nFyR5TN+yof3O5lRVTvOY6HW0XwfsA2wPfA94/GZt3gx8vM2vBL7Q5vcHdm/zTwQ29K1zIbBijLktA66YZbsXA88AApwHvHCUuW3W5knAdWPYb8uAJwOnAS/ri+8CXN9ed27zOy/EftvKvPYD9m3zuwMbgUe096f2tx31PmvLfjHLds8EVrb5jwNvGmVem32vPwEeMuJ9dlDfZ76Je/77HNrvrGvyyGL+fjsESVX9P2B6CJJ+hwOfafNfAp6TJFX13ar6UYtfCTw4yQ6TkNtsG0yyG7BTVV1UvV/macARY8ztyLbuQurMrapurKq1wG82W/cFwPlV9ZOq+r/A+cDBC7Tf5p1XVf2gqq5t8z8CbgGmtvDzh5LbbNp3/Wx63z30fgtHjCmvlwHnVdWvtvDztza3f+z7zIvo3VsGw/2dzcliMX8zDUGyx2xtquou4HbgkZu1eSlwWVXd0Rf7dDvE/c/zPJTc2tz2TvLdJN9M8qy+9us7tjmK3Ka9Ajh9s9go9tuWrrsQ+21r8vqtJAfQ+0v2ur7wCe1Ux4fn+QfL1ub2oCSrk1w0faqH3nf90/bdz2ebC5HXtJXc93c26n32OnpHCnOtu1D/fc7KYjFGSZ4A/CXwhr7wH1fVk4BntelVI05rI7BXVe0PHAd8PslOI85hTkmeDvyqqq7oC497v0209pfnZ4HXVNX0X9LvBH4X+D16pzXeMYbUHlO9ISz+HfCRJI8dQw4zavvsSfTu75o20n2W5JXACuCvh/k5g7BYzN8gQ5D8tk2S7YCHA7e190uBs4Gjquq3f+lV1Yb2+nPg8/QOWUeWW1XdUVW3tRwupfdX6H6t/dK+9ec75MpW7bfmPn/tjXC/bem6C7Hftmq4m1bsvwa8u6oumo5X1cbquQP4NKPfZ/3f2/X0+p32p/ddP6J991u8zYXIq3k5cHZV3dmX78j2WZLnAu8GDus78zDM39ncFrIDZDFN9O5+vx7Ym3s6qZ6wWZtjuXdH7Zlt/hGt/Utm2Oaubf6B9M7ZvnHEuU0BS9r8Pu0Ht0vN3IF2yChza+8f0HLaZxz7ra/tqdy3g/sGep2OO7f5BdlvW5nX9sAFwFtnaLtbew3wEeDEEe+znYEd2vyuwLW0jl7gi9y7g/vNo8qrL34RcNA49hm9onkd7eKEUfzOOvNeyI0ttgk4BPhB+1Lf3WIfoPeXAMCD2o9+Xfsi92nx9wC/BNb0TY8CHgpcCqyl1/H9Udr/uEeY20vbZ68BLgP+Td82VwBXtG1+jDYCwKhya8sOBC7abHuj3G+/R+988C/p/QV8Zd+6r205r6N3umfB9tt88wJeCdy52W9teVv2P4HLW25/BzxslPsM+IP2+d9rr6/r2+Y+7btf134LO4z4u1xG74+SB2y2zVHts28AN/d9Z+eM4nc21+RwH5KkTvZZSJI6WSwkSZ0sFpKkThYLSVIni4UkqZPFQpLUyWIhzVOSpyT5VhtK+jdJKskHxp2XNAzeZyHNQ5IH0btZ6qiqujjJn9O7mfDt5X9Uuh/yyEKan+fSGy344vZ+Lb1hFyrJJ4f1oQv+QBtpQBYLaX6eSG/Yh2lPBS5L8mDgXyd5X5IzkjwwyX9J8pEkJwEkmUry6SRLk5yS5D8kOagtOyXJw9r6H03yZ0l+J72n3b0D+J2R/0slegNaSdpyt9F7QA9J9gNeQm+so/2BL1bVR5L8d+DtwIOBn9IbOI6q2pTkn4EP0ntWwXLgKUnuBi6hN8rudm2dZ9IbI+n0qjppRP826T4sFtL8nA4cluQK4FbgyKq6rT1gaG1r8xB6BeLY6nu4VZKH0Rso766q+kXbxpH0Brb798Aq4E/pjQC8J71i8uWR/KukWdjBLS2gJJ+gPY+B3uMyd6T3YJ+b6I1Y+g16xeD99J6XcElVXZjkcuC9VfX3Sd4G7ETviXHX0CsWr697HlokjZzFQpLUyQ5uSVIni4UkqZPFQpLUyWIhSepksZAkdbJYSJI6WSwkSZ0sFpKkThYLSVKn/w/Mcv0r9yyyzQAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "plt.hist(bexvar_distribution, bins=20)\n", "plt.ylabel(\"# of samples\")\n", "plt.xlabel(r\"$\\sigma_{bexvar}$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If the light curve is intrinsically variable, then the posterior distribution of bexvar should exclude low values. Users can compute \n", "the lower 10% quantile of the posterior, and use it as a variability indicator (see [Buchner et al. (2021)](https://arxiv.org/abs/2106.14529))." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The method uses fractional exposers (`frac_exp`) in each bin to compute the count rates (i.e.$~\\scriptstyle{R_i = C_i/(\\Delta{t_i}\\times f_i)}$). In its current form it only considers time bins with `frac_exp` < 1. The `bg_ratio` parameter is used to scale the `bg_counts` to estimate counts in source region. The `bg_count`, `bg_ratio` and `frac_exp` are optional parameters, if they are not provided, the method defines default values for them as described in documentation. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us see an example to get bexvar distribution without these optional parameters." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "preparing time bin posteriors...\n", "running bexvar...\n", "[ultranest] Sampling 400 live points from prior ...\n", "[ultranest] Explored until L=-4e+01 [-36.8486..-36.8486]*| it/evals=3615/5101 eff=76.8985% N=400 \n", "[ultranest] Likelihood function evaluations: 5125\n", "[ultranest] logZ = -41.34 +- 0.09729\n", "[ultranest] Effective samples strategy satisfied (ESS = 1692.4, need >400)\n", "[ultranest] Posterior uncertainty strategy is satisfied (KL: 0.46+-0.06 nat, need <0.50 nat)\n", "[ultranest] Evidency uncertainty strategy is satisfied (dlogz=0.10, need <0.5)\n", "[ultranest] logZ error budget: single: 0.09 bs:0.10 tail:0.01 total:0.10 required:<0.50\n", "[ultranest] done iterating.\n", "\n", "logZ = -41.331 +- 0.174\n", " single instance: logZ = -41.331 +- 0.092\n", " bootstrapped : logZ = -41.335 +- 0.174\n", " tail : logZ = +- 0.010\n", "insert order U test : converged: True correlation: inf iterations\n", "\n", " logmean : -0.517│ ▁ ▁ ▁▁▁▁▁▁▁▁▁▁▂▂▃▄▆▇▇▆▅▄▃▂▁▁▁▁▁▁▁▁▁ │0.383 0.020 +- 0.081\n", " logsigma : 0.029 │ ▁▂▆▇▇▅▃▂▂▁▁▁▁▁▁▁▁ ▁▁▁▁▁ ▁ ▁ │1.236 0.213 +- 0.074\n", "\n", "running bexvar... done\n" ] } ], "source": [ "import numpy as np\n", "\n", "time = np.arange(0,8)*100\n", "counts= np.array([106, 87, 115, 148, 43, 129, 204, 87])\n", "time_del = np.ones(np.size(time))*100\n", "\n", "bexvar_distribution = bexvar.bexvar(time=time, src_counts=counts, time_del=time_del)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Bexvar: Theoretical background" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "This section provides a theoretical understanding of Bayesian excess variance (bexvar). This is an optional read.\n", "\n", "Given a lightcurve data ${\\scriptstyle 𝐷 = (𝑆_1,𝐵_1,~…~,𝑆_𝑁,𝐵_𝑁)}$\n", " where ($\\scriptstyle{S_i}$) denotes counts obtained from source region and ($\\scriptstyle{B_i}$) denotes counts obtained from background extraction region in $\\scriptstyle{i^{th}}$ time bin.\n", "If it is assumed that the counts $\\scriptstyle{𝑆_𝑖}$ and $\\scriptstyle{𝐵_𝑖}$ can be expressed as\n", "Poisson processes. \n", "$$ \\scriptstyle {𝑆_𝑖 ~ \\sim ~ 𝑃𝑜𝑖𝑠𝑠𝑜𝑛((( 𝑅_𝑆(𝑡_𝑖) ~+~ 𝑅_𝐵(𝑡_𝑖) \\times 𝑟)~×~𝑓_𝑖~\\times~Δ𝑡)}$$\n", "$$ \\scriptstyle {𝐵_𝑖 ~ \\sim ~ 𝑃𝑜𝑖𝑠𝑠𝑜𝑛(𝑅_𝐵(𝑡_𝑖) × 𝑓_𝑖 × Δ𝑡)}$$\n", "\n", "Here, $\\scriptstyle{𝑅_𝑆(𝑡_𝑖)}$ is source count rate and $\\scriptstyle{𝑅_B(𝑡_𝑖)}$ is background count rate in $\\scriptstyle{i^{th}}$ time bin. \n", "It is further assumed that $\\scriptstyle{𝑅_𝑆(𝑡_𝑖)}$ is distributed according to a log normal distribution, with some unknown parameters (i.e., $\\scriptstyle{log(\\bar{𝑅_{S}})}$, and $\\scriptstyle{\\sigma_{bexvar}}$).\n", "$$\\scriptstyle{log(𝑅_𝑆(𝑡_𝑖))~\\sim~𝑁𝑜𝑟𝑚𝑎𝑙(log(\\bar{𝑅_𝑆}),~ \\sigma_{𝑏𝑒𝑥𝑣𝑎𝑟})} $$\n", "\n", "This $\\sigma_{𝑏𝑒𝑥𝑣𝑎𝑟}$ provides intrinsic variability on log-count rate and it is defined as Bayesian excess variance (bexvar). The posterior distribution of $\\sigma_{𝑏𝑒𝑥𝑣𝑎𝑟}$ can be used to identify intrinsically variable object.\n", "\n", "The bexvar() method in Stingray returns posterior samples of $\\scriptstyle{\\sigma_{𝑏𝑒𝑥𝑣𝑎𝑟}}$ given a light curve data.\n", "The samples are generated following the same prescription given in [Buchner et al. (2021)](https://arxiv.org/abs/2106.14529). The method uses flat, uninformative priors on $\\scriptstyle{log(\\bar{𝑅_𝑆})}$ and $\\scriptstyle{log(\\sigma_{𝑏𝑒𝑥𝑣𝑎𝑟})}$ and obtains the posterior samples using nested sampling Monte Carlo algorithm MLFriends (Buchner [2016](https://link.springer.com/article/10.1007/s11222-014-9512-y),\n", "[2019](https://arxiv.org/abs/1707.04476)) implemented in the [UltraNest](https://johannesbuchner.github.io/UltraNest/) Python package (Buchner [2021](https://arxiv.org/abs/2101.09604)).\n", "\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" }, "vscode": { "interpreter": { "hash": "f6246b25e200e4c5124e3e61789ac81350562f0761bbcf92ad9e48654207659c" } } }, "nbformat": 4, "nbformat_minor": 4 }